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Abstract 

We suggest a method to create quantum turbulence (QT) in a trapped atomic Bose-Einstein 
condensate (BEC). By replacing in the upper half part of our box the wave function, ^, with its 
complex conjugate, ^* , new negative vortices are introduced into the system. The simulations 
are performed by solving the two-dimensional Gross-Pitaevskii equation (2D GPE). We study the 
successive dynamics of the wave function by monitoring the evolution of density and phase profile. 

PACS numbers: 03.75.Kk, 03.75.Lm, 47.27.tb 
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I. INTRODUCTION 



The aim of this paper is to explore the physics of a turbulent vortex system. Turbulence 
is studied both in classical and in quantum fields. 

A dynamical and statistical description of the classical turbulent flow has been provided 
earlier by Kolmogorov [1']. The turbulent flow is characterized by random and chaotic 
three-dimensional vorticity and irregularity. Rapid mixing of momentum, heat and mass 
is manifested as diffusivity, which is a key characteristic of the turbulent flows. Under the 
influence of viscosity in the turbulent motion, the kinetic energy is dissipated into heat. 
Turbulence happens with high Reynolds numbers and is generally anisotropic. Simulations 
using the Gross-Pitaevskii equation (GPE), modelling qualities of both classical and quan- 
tum turbulence, have been well studied 

Key properties of superfluid vortex lines were discovered by Onsager, 5| and later de- 
veloped by Feynman, 6[. In superfluids like ^He or ^He, QT was studied with quantized 
vortices characterized by topological defects Q]- Conflgurations of quantized vortices were 
investi gate d in superfluid ^He js], which can be g rouped into two types [o^: ordered vortex 



arrays 



10| and disordered vortex tangles 



Disordered vortex tangles appear, when 



lelium is made turbulent under thermal counterflow velocity or using^ grids or propellers 
I3I . for example by moving a grid at T = O.OlTx where, Tx = 2.17 [14] or by a vibrating- 



wire resonator in B-phase superfluid ^He (^He-B) [15 



Large scale turbulence of quantized 
8|. The disadvantage of turbulence in 



vortices was studied in superfluid ^He-B [TJ and ^He 
BEG is the small system size and the relatively low number of vortices. On the other hand, 
the advantage is the relatively good visualization in detail of the individual vortices. 

In magnetically or optically trapped rotating atomic Bose-Einstein condensates (BEGs), 
vortex lattices with quantized vortices along the rotation axis are formed. The density and 
phase of BEGs can be directly observed. The rotational motion is sustained by quantized 
vortices and is quantized in units of /t = q{h/m), where q is integer, m is the particle mass 
and h is the Planck constant. For g > 1 the system is unstable . 

The circulation around each superfluid vortex fllament is flxed by the condition: 
Vsdr = K where, C is a circular path around the axis of the vortex and is the local 
velocity flow. 



Quantized vortices are characterized by singularities in the phase and by associated holes 
in the density. Vortex lattices were first realized experimentally by the groups of Madison 
[l6 | and Abo-Shaeer [l^]. The corresponding simulations were made for example, by the 
group of Tsubota [isl - 

Similarities between QT and classical turbulence (CT) were observed experimentally 
|19l . 120| . A realizable study of QT was presented by using numerical simulations of the GPE 



in trapped BEG, and by combining rotations around two axes [2l|]. The kinetic energy, 
i^kin, can be divided into a compressible part, Ef^^, due to the sound waves, and into an 
incompressible part, E^^^, due to the vortices [3]. It was shown that, the spectrum of 
the incompressible kinetic energy obeys the Kolmogorov law and the energy flux becomes 
constant value. On the other hand, the method we present here is a pure mathematical 
model, which can give ideas and inspirations to the experimental investigations. With this 
method a superfiuid analogy of classical 2D turbulence was created in a 2D system. 

The wave function or the order parameter of the solution of GPE is a complex number, 
ip, corresponding to the amplitude of the particle to have a given position, r at any given 
time, t. The condensate wave function, ip in a complex plane or Argand diagram is observed 
as a positive vector with the real part. Re ip and the imaginary part, Im ip. The real and 
imaginary parts are regarded as independent quantities. 

The condensate density is defined as n{r, t) = \ip{r, t)p, while the phase of the condensate 
as S" = tan~^ [Im '?/'(r, t)/Re '?/'(r, t)]. Using the Madelung transformation, that is to say 
the polar form of in terms of the superfiuid density and macroscopic phase, we obtain 
ip{r,t) = ^/\ip{r, t)|2e*'^*^'"'*) = ^/7^Jr~t)e^^'''^'*\ The superfiuid local velocity flow is given by 
v,(r,t) = {h/m)VSir,t). 

The complex conjugate of a complex number is given by changing the sign of the imag- 
inary part. The complex conjugate of '?/'(r, t) is ip* = lie ip — i Im ip and similarly 
ip*{r,t) = y^[^^^(rvOF^~*'^^'''*'' — V '^('^5 t)e~*'^^'^'*''. In the complex plane, ip* is symmetric 
about the real axis and ip + ip* and ip ■ ip* are real numbers. If a complex number supplies 
a solution to a problem, so its conjugate does too. 
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II. THEORY: THE GROSS-PITAEVSKII EQUATION 



To explore the dynamics of superfluid vortices at nonzero temperatures and the density 
and phase profile of a rotating condensate, we solve numerically the following 2D GPE for 
the condensate macroscopic wave function, ip{r,t): 



(9^ 

dt 



2m 



Vi + Krap + g2DN\^\^ - fl - QL, 



(1) 



where, m is the atomic mass, a is the s-wave scattering length, Q is the angular frequency of 
rotation about the z-axis and Lz = ih{xdy — ydx) is the angular momentum operator. The 
rotation frequency, Q is large enough, so several vortices are present in the field, ip forming 
a regular array. The phenomenological damping parameter is 7 = 0.01, which models the 
interaction of the condensate with the surrounding thermal cloud 23|, |2J]. In 2D due to the 



enough strong confinement along the z-axis the coupling constant becomes 



22 



g2D = V327ThQ 



L 



(2) 



By using the following formulas, = h/2mQ and = h/muz and by multiplying the 
numerator and the denominator by Uz, we obtain finally: 



g2D = 2V2nhalzUJz, (3) 

where, Q is the angular velocity along the z — axis, I is the magnetic length, lz is the 
characteristic length of the z-axis oscillator and uJz is the characteristic frequency of the 
trapping z potential. In dimensional form the chemical potential, /x is defined as: 

In imaginary time, excitations are damped and both ip{r,t) and /i converge to a stationary 
solution supplying exact initial conditions for time-dependent solutions. In 2D the harmonic 
trapping potential is defined by: 

Krap(a;, y) = ^mul {x^ + y^) , (5) 
here, the radial trap frequency is uj±. The ground state properties of disk-shaped BECs 
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{ujz ^ uj±) are characterized by only one parameter 



25| 



Na 



(6) 



here, is the number of particles, a is the s-wave scattering length, is the axial oscillator 
length, and A = u!z/uj± ^ 1, is the trap aspect ratio. With the help of A, Iz and l± = 



yjhjwjjj^ (which is the radial oscillator length) expressions, we obtain the final form of K^. 



K2 = Nall/li. (7) 

The perturbative regime (quasi-2D regime) corresponds to K2 <^ 1, while for the Thomas- 
Fermi (TF) regime we have the condition, K2 ^ 1. In the perturbative regime, the axial 
wave function coincides essentially with the ground state (Gaussian) wave function of the 
corresponding (z) harmonic oscillator. On the other hand, in the TF regime the axial wave 
function is essentially a TF wave function. 

Thus, only when K2 ^ 1 one can assume that, the wave function along the 2;-axis is the 
ground state of the harmonic potential. This only occurs when, and hence the nonlinear 
mean-field interaction term in the equation of motion is small enough or A is large enough. 
There exists a direct relation between K2 and the parameter, 0/2^1-2(0) where, n ofr 1 , t) is 
the local condensate density per unit area characterizing the radial configuration 261]. 



n2{T^,t)=N j dz\^ir^,z,t)\^. (8) 

In fact, for K2 <^ 1 ^ 0/2^0 <^ 1 and for K2 ^ 1 ^ a/^no ^ 1. The (axial) local chemical 
potential, ftz = l^'zlhjJz [26] is used in [27] to obtain an effective 2D equation of motion for 
disk-shaped condensates. In the quasi-2D mean- field regime, K2 ^ 1 ^ alzn2 ^ 1. This 
effective equation reduces to Eq. [H where A^|\l/(0)p = ^2(0) ~ n2. One cannot know the 
value of n2 in advance, therefore in practice K2 is the relevant parameter, which is a global 
parameter. In the TF regime, ^ 1 ^ alz'n2 ^ 1- The effective 2D equation |27| reduces 
to: 



^^^^ 



/ 37r ^ 



(9) 
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In the above equation, we are assuming \1/ to be normalized to unity. Both Eq. [T] and Eq. [9] 
are the correct hmits of the underlying three-dimensional GPE. 



It is also necessary that, the typical time scale of the radial motion, (At) to be much 
larger (i.e. slower) than the time scale of axial motion (which is of the order of uj~^). This 
is necessary, in order for the adiabatic approximation to be valid, and as a result radial and 



axial motions to be separable ! 



27|. 



For stationary problems, ~ oo, and for collective oscillations, A^ ~ and thus the 
condition, A = uJz/uj± ^ 1 guarantees the fulfillment of the adiabatic approximation. Con- 
nect with Eq. [H there is no problem, since A^ ~ {ujzalzn2)'^ ^ ^7^? which is a consequence 
of the fact that, in this case 0/2^2 ^ 1. 

The total energy, -Etot? can be identified with the sum of three energies: 

Etot = E]^m + Eint + -Strap, (10) 

where, the kinetic energy, -Ekin, the internal energy, Ei^^t, and the trap energy, -Etrap are given 
respectively by: 

EUt) = J |^(Vp(x,t)v(x,t))'d\ (11) 
Eintit) = i y g{p{^,t)fd\ (12) 

^trap(t) = j pix,t)VUx)dh. (13) 

In the kinetic energy expression, the main contribution to he energy comes from the phase 
gradients. The contribution from the density gradients is neglected due to its small value. 

It is useful to scale the GPE in dimensionless units. We use harmonic-oscillator units 
(h.o.u.) in the whole paper where the units of time, length and energy are: ujj_ , 

y^h/muj_i and huj± respectively, so that: 



7) 



dip 
'dt 



-Vi + Vtrap + C|^|2-/i-^]L, 



(14) 



where, 
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^trap = ^(a;' + y'), (15) 

and the effective coupling constant, which is an important characteristic parameter of the 
2D system becomes: 

where, g^D is expressed by Eq. |2] and by Eq. [31 Using the expression from Eq. [3] together 
with the expression of /_l = a/ ^/ muo±^ and then by substituting muOz/h by l//^, we obtain 
first: 

C = 2v^^, (17) 

and later by substituting Na/l^X^ by -^'2, we obtain in 2D the final form of the effective 
coupling constant: 



C = 2V27iX^K2. (18) 

Here, N represents the number of atoms per unit length along z. Throughout this paper we 
use, C = 1400 in these calculations. In order for Eq. [T3]to be valid, we need K2 <^ 1. Say 
for example, K2 = 0.1: 

which, implies that: 

» (20) 



and 



u^l » ^^i- (21) 
Vovr 



We would like to present an example for the experimental realization. For a ^''Rb con- 
densate, the s-wave scattering length is a = 5.82 nm, and the axial oscillator length is: 
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/ ^ 10.78 

= \ = ^^/^"^- (22) 



Substituting these values for a and Iz in Eq. [TTl we obtain: 



C = 2.7066 X 10-^vW27rA^, (23) 

and using, C = 1400 one finds: 

517254.03 , 

N= (24) 

Considering typical values, as for instance, ujz/2n = 600 Hz and LJ±/2n = 10 Hz, one has a 
condensate with = 21117 ~ 21100 particles and: 

A = uju^ = 60, (25) 

and 

C 

K2 = ^ = 0.0776 ~ 0.08. (26) 
Moreover, one can easily estimate (rather accurately) the radius, R and the chemical poten- 



tial, /i of the ground state of this condensate |25|: 



/i 



^ + (2^/2/^X2) = 0.852, (27) 



fkUz 2 
and 

^ = VX ^'^ = 6.498. (28) 

Using that, = ^/Xl^ = y/oo 10.78/7600 = 3.41 jjm, we find R = 22.15 /im. On the other 
hand, since fi = 0.852 huj^ < f^z, we conclude that, the condensate is indeed in its axial 
ground state (quasi-2D regime). However, 

^ A-^ = 51.12, (29) 



hu!± huz 

^ yU ^ huj±, which means, that with respect to its radial motion, the condensate is in a 
Thomas- Fermi regime (many radial modes excited). This is also confirmed by the fact that, 
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We consider a strongly interacting 2D BEC in a harmonic trap rotating at the angular 
frequency, Q. The numerical calculations are performed with the semi-implicit, Crank- 
Nicholson method in a square box of size +/— 7/8 (h.o.u.). 

III. RESULTS 

First, we create a non-rotating condensate at t = 0, and later Q is set to 0.85/0.8 to 
create a stable lattice of 22/20 vortices in a rotating frame. Physical quantities of the BECs 
like density and phase can be clearly observed. Arrays of vortices are shown by the density 
profile of the condensate for Q = 0.85 at t = 199, see Fig. [T] (up) or by the density and phase 
profile of the condensate for Q = 0.8 at t = 199, see Fig. [2l Rotation of a superfluid takes 
place via Abrikosov lattice of quantised vortices, where the rotational velocity profile mimics 
solid body rotation. Away from vortex cores the superfluid is irrotational. For vortex lattice 
with vortices, the circulation becomes, F = J ■ (il = Nk, and the vortex density depends 
only on Q, as we can see from its expression, riy = N/A = 2Q/n. Here, fl is the angular 
velocity of the trapping potential and A^ is the number of vortices. The circulation of the 
fluid is quantised in units of k. 

Disordered vortex arrays are produced with the help of the following method: \1/ is 
instantaneously changed to \E'* in the upper half part of the box {y > 0; the center of the 
box is at [x,y] = [0,0]). That means that, the condensate is divided into two equal parts 
with opposite rotation. The upper part rotates anti clockwise and the bottom part rotates 
clockwise. These anti-vortices in the upper part of the condensate change their sense of 
rotation and are visualized by Fig. [1] (down). 

This method is a shock for the condensate. First, the largest part of the upper side of 
the condensate becomes deformed and after some shaking movements, see Fig. |3l the system 
recovers its circle-like form again with turbulence, shown by Fig. |H So, the method can thus 

)e used to create a 2D turbulence in the same way as in the classical Onsager vortex gas 

29|. 

After using of this method, we distinguish two different cases: [fl = 0.8) and ($7 = 0). 
We analyze the total number, A^ of vortices, as well as the number A^+, of positive and A^~, 
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of negative vortices (see Fig. [6] and Fig. [7]) under certain conditions, for example keeping Q 
constant or not. For Q = 0.8, after a certain transition period, the system enters the state 
which is the same as the original regular vortex lattice. 

IV. DISCUSSION 

First, we create vortex lattice in a rotating trapped BEC. Then, we change the wave 
function with its complex conjugate by instantaneously reversing the direction of rotation 
of the atoms in the top half part of the box. As a consequence, in the upper part of the 
trap, we transform vortices into their anti-vortices rotating in the opposite direction. 

Thus, we suggest a method to create a system of positive and negative vortices in a disk- 
shaped condensate, which is divided into two equal parts rotating in opposite directions, 
as shown in Fig. [1] (down). To prepare the right conditions for turbulence, at t = 200, the 
rotation frequency, Q was set to 0. As a result, the kinetic energy, E'km and the 2;-component 
of the angular momentum, tend to zero, (see. Fig. [5l). 

In that case, the negative vortices remain in the condensate, as seen in Fig. [3 We find 
that, after applying the recent method, the new anti- vortices interact with the existing 
positive vortices and contribute to the formation of turbulence, see Fig. HI 

In two separate cases, we track the number of vortices and anti-vortices by maintaining 
or discontinuing the trap rotation. If we continue to use Q = 0.8 or 0.85 after the suggested 
method, the anti-vortices move out to the edge of the condensate, and new vortices come in 
until the vortices settle into an ordered vortex array. The corresponding total number, N, 
positive, , and negative, N~ , numbers of vortices are shown in Fig. El when the system 
will relax to the original vortex array. 

On the contrary, in the second case (for i7 = 0) both vortices and anti- vortices coexist and 
interact with each other to create turbulence. In that case, more negative vortices and more 
possibilities for annihilations appear, see the suitable plots for A^, A^"*" and A^~ in Fig. [71 
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V. CONCLUSION 



It is instructive to compare this method (call A) to create turbulence with complex 



conjugation of the order parameter in a semi-plan with a different one (call B) 30|, when 
the phase was imprinted in upper left quadrant (x < and y > 0) and bottom right quadrant 
(x > and y < 0) of the box, or with a method to produce positive and negative vortices 
by the phase imprinting method (in upper two quadrants), when soliton-like perturbation 



due to the snake- instability, 3l| decays into vortex - anti vortex pairs 32] • In that case, the 



total number of negative (anti)vortices is not enough to create turbulence. 

In case B, the direction of vorticity of quantized vortex-line (in 2D vortex-point) or the 
sign of circulation does not change. On the other hand, in case A the direction of the 
condensate does flip. Thus, in B, we change the sign of the square root of the density in 



the Madelung transformatio: ip{r,t) = —^y\ip(r, t) ^e*'^*^'"'*) = — A/ri(r, t)e*'^*^'''*\ Whereas, in 
A, we change the sign of the phase: ip{r,t) = -^/|^^(ivt)pe~*'^*^'''*) = ^/n{r, t)e~*'^*^'"'*^. 

Another important difference between these two methods is that, in the case B soliton- 
like perturbations appear along the x- and y-axis. The solitary wave with local density 
minimum and a sharp phase gradient of the wave function at the position of the minimum, 
embedded in a 2D geometry leads to dominant decay mechanism. Thus, the soliton-like 
perturbations bend and decay into a more stable vortex-anti vortex pairs accompanied by 
sound waves. In our case A, we do not observe soliton-like perturbations! Furthermore, 
after applying these two methods {A and B), the shape of the condensate is different, due 
to the fact that, in A we do not have solitary waves. 

Application of these two methods {A and B) lead to different physical behaviour. For 
example, the number of total, positive and negative vortices as a function of time show 
different behaviour in A and in B. With the method A more negative vortices appear, 
which contribute to more annihilations and different total number of vortices as a function 
of time. Compare, for example these numbers of vortices, for the case when Q was set to 
after the application of these methods. Our observation is that, the variation of different 
numbers of vortices with time vary more slowly in the case A. 

We conclude that, in case A the complex conjugation method reverses the rotation of the 
vortices in the two halves of the box. Now, the system contains roughly the same number 
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of vortices with different sign and witfi different energies due to the corresponding different 
velocity fields. Thus, the imitation of the sohd body rotation becomes damaged (stop of 
solid body rotation and far field almost cancels). 

The phase across the whole condensate follows the sign of the vortices, and becomes 
different in the upper and lower parts of the box. At the centre of the box, the phase is zero 
and by that, no phase kinks happen upon complex conjugation. Thus, no soliton like per- 
turbations form. Very different energies appear and contribute to difficulty in experimental 
utilization. We do not know a way to realize this method experimentally. 

Our study about a more global properties of the phase across the whole condensate 
shows that, this new method of complex conjugation change the rotation of the vortices in 
the upper half part of the box. On the other hand, with the phase imprinting method the 
system try to smooth out the change caused by generation of a discontinuity in the phase 
in the form of solitary and sound waves. 
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FIG. 1: For = 0.85, the density profile of the condensate at t = 199 (left) and t = 200 (right), 
after applying the method described in this paper. The vortices rotate in opposite directions in the 
upper half and in the bottom half part of the condensate. (The squared shape of the condensate 
depends on the small size and the boundary effect of the box.) 
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FIG. 3: The density (left) and phase (right) profiles of the condensate, corresponding to FIG. [2] at 
t = 202.1. 
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FIG. 5: Different energies and for the presented method, as a function of time, t (from bottom to 
top: kinetic energy (green), z-component of the angular momentum (red), internal energy (blue), 
trap energy (purple) and total energy (cyan)). To keep the anti- vortices inside of the condensate 
after t = 200, was set to 0. 
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FIG. 6: The number of positive vortices (left), and (right) the total number of vortices (A'^, 
represented by triangles) along with the number of negative vortices (A^~, represented by circles) 
for the presented method. Here, was kept at 0.8 after t = 200. 
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FIG. 7: The total number of vortices (A^, represented by triangles), the number of positive vortices 
(A^+, represented by circles) and the number of negative vortices (A^~, represented, by filled circles), 
corresponding to FIG. [6] for the case of O = from t = 200 (when the method was applied). 
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